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Abstract 

We present the TRIQS library, a Toolbox for Research on Interacting Quantum Sys¬ 
tems. It is an open-source, computational physics library providing a framework for 
the quick development of applications in the field of many-body quantum physics, and 
in particular, strongly-correlated electronic systems. It supplies components to develop 
codes in a modern, concise and efficient way: e.g. Green’s function containers, a generic 
Monte Carlo class, and simple interfaces to HDF5. TRIQS is a C++/Python library 
that can be used from either language. It is distributed under the GNU General Public 
License (GPLv3). State-of-the-art applications based on the library, such as modern 
quantum many-body solvers and interfaces between density-functional-theory codes and 
dynamical mean-field theory (DMFT) codes are distributed along with it. 
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Nature of problem: 

Need for a modern programming framework to quickly write simple, efficient and higher-level 
code applicable to the studies of strongly-correlated electron systems. 

Solution method: 

We present a C++/Python open-source computational library that provides high-level abstrac¬ 
tions for common objects and various tools in the field of quantum many-body physics, thus 
forming a framework for developing applications. Running time: Tests take less than a minute. 
Otherwise highly problem dependent (from minutes to several days). 


1. Introduction 

In this paper, we present the 1.2 release of the TRIQS project (Toolbox for Research 
in Interacting Quantum Systems), a free software library written in Python and C++ for 
the implementation of algorithms in quantum many-body physics. TRIQS is distributed 
under the GNU General Public License (GPLv3). 

Strongly-correlated quantum systems are a central challenge for theoretical condensed 
matter physics with a wide range of remarkable phenomena such as metal-insulator tran¬ 
sitions, high-temperature superconductivity and magnetism. In the last two decades, 
tremendous progress has been made in the field of algorithms for the quantum many- 
body problem, both in refining existing techniques and in developing new systematic 
approximations and algorithms. Methods to address the quantum many-body problem 
include dynamical mean-held theory (DMFT) [1, 2] and its cluster [3] or diagrammatic 
extensions [4, 5] or the density matrix renormalization group (DMRG) [6] . DMFT meth¬ 
ods can also now be combined with more traditional electronic structure methods such as 
density functional theory (DFT) leading to ab initio realistic computational techniques 
for strongly-correlated materials [2]. Several collaborative software development efforts 
have made some of these theoretical developments largely accessible, e.g. Refs. [7, 8, 9]. 

The purpose of the TRIQS project is to provide a modern framework of basic building 
blocks in C++ and Python. This is needed for the rapid implementation of a broad spec¬ 
trum of methods. Applications range from simple interactive phenomenological analysis 
in Python to high-performance quantum impurity solvers in C++. At this stage, TRIQS 
focuses primarily on, but is not limited to, solid-state physics computations, diagram¬ 
matic approximations and methods of the DMFT family (DMFT, clusters and underlying 
quantum impurity solvers). 

A particular emphasis is placed on the documentation, in particular in providing short 
code examples that can be reused immediately (in Python and C++). Full documentation 
of the project is available online: http://ipht.cea.fr/triqs. 

Several applications have already been built with the TRIQS library, and some of 
them are publicly distributed. Let us mention a state-of-the-art implementation of the 
hybridisation-expansion quantum impurity solver CTHYB (http: //ipht. cea. fr/triqs/ 
applications/cthyb/) and the DFT TOOLS project which provides an interface be¬ 
tween DMFT and DFT packages such as WIEN2k for realistic computations for strongly- 
correlated materials (http://ipht.cea.fr/triqs/applications/dft_tools/). Since 
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these applications are not part of the library itself and involve a different set of authors, 
they will be presented in separate publications. However, they are distributed along 
with the TRIQS library under GPL license and are available for download on GitHub 
(https://github.com/triqs) . 

The TRIQS project uses professional code development methods to achieve the best 
possible quality for the library and the applications, including: i) version control using 
git; ii) systematic code review by the main TRIQS developers; in) test-driven develop¬ 
ment: features of the library are first designed with a series of test cases. When the 
implementation is completed, they become the non-regression tests executed during the 
installation process. 

This paper is organised as follows: we start in Sec. 2 with the main motivations for 
the project. In Sec. 3, we outline the structure of the TRIQS project. Sec. 4 summarizes 
our citation policy. In Sec. 5, we discuss the prerequisites to efficient usage of TRIQS 
and Sec. 6 describes the portability of the library. In Sec. 7, we provide two illustrating 
examples to give a flavour of the possibilities offered by the library: we show that TRIQS 
makes it possible to write a DMFT self-consistency loop in one page of Python, and, in 
another example, how equations can be coded efficiently in C++. In Sec. 8 we review 
the most important library components. In Appendix A we present the implementa¬ 
tion of a fully working, MPI-parallelized, modern continuous-time quantum Monte Carlo 
algorithm (the so-called CT-INT algorithm [10, 11]) in about 200 lines. This example 
illustrates how TRIQS allows one to design a complex, yet short, readable and extensible 
code. 

2. Motivations 

The implementation of modern algorithms for quantum many-body systems raises 
several practical challenges. 

Complexity: Theoretical methods and algorithms are becoming increasingly complex 
(e.g. quantum Monte Carlo [10, 12, 13, 11] and dual boson [14] methods). They are hence 
more difficult to implement, debug and maintain. This is especially true for realistic 
computations with methods of the DMFT family, where one has to handle not only the 
complexity of the many-body problem but also the various aspects (orbitals, lattices) of 
real materials, which usually requires a well-organised team effort. 

Versatility/Agility: Algorithms are changing and improving rapidly, sometimes by 
orders of magnitude for some problems [15]. This can lead to a possibly quick obsoles¬ 
cence of a given implementation. To address new physics problems requires regular and 
substantial modifications of existing implementations. Moreover, there are numerous 
possibilities for new algorithms which need to be tested quickly. 

Performance: Modern algorithms are still quite demanding on resources, e.g. quantum 
Monte Carlo methods. Hence, the performance of the codes is critical and a simple 
implementation in a high-level language is usually not sufficient in practice. 

Reproducibility: The central role and the growing complexity of the algorithms in our 
field reinforces the need for reproducibility, which is central to any scientific endeavour. 
Therefore, the results obtained by a numerical computation should in principle be pub¬ 
lished systematically along with the code that produced them, in order to allow others 
to reproduce, falsify or improve on them. This requires codes to be readable (i.e. written 
to be read by other people than their author) and relatively quick to produce. 
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To address these challenges, one needs readable, clear and simple implementations 
with reusable components provided through high-level abstractions. We emphasize that 
this is not in contradiction with the requirement of high performance. The combination 
of a higher level of abstraction with high performance is achieved using modern program¬ 
ming techniques (e.g. generic programming). The purpose of the TRIQS project is to find 
and efficiently implement the relevant abstractions, basic components and algorithms for 
our domain. 


3. Structure 

The TRIQS framework is depicted in Fig. 1: the core library is at the root (bottom) 
consisting of basic building blocks, which are used in a series of applications (top). The 
applications can be in pure python (e.g. DFT tools), in C++ with a Python interface 
(e.g. CTHYB), or even in pure C++. The subject of this paper is the core library. 

Figure 1: Structure of the TRIQS project 
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The components of the TRIQS library can be used both in Python and in C++: C++ 
brings the performance needed for applications where speed is critical (like many-body 
solvers) and the type safety of a compiled language. On the other hand, Python is typi¬ 
cally used as a higher-level interface for data analysis, investigation of phenomenological 
approaches, and tasks related to reproducibility. Most objects of the library are written 

4 


































in C++ and exposed to Python using a specially designed tool described in Sec. 8.8. As 
a result, TRIQS can be used together with all the modern scientific tools of the Python 
community, in particular with IPython notebooks [16] which are recommended for an 
optimal interactive usage of the library in Python. 

4. Citation policy 

We kindly request that the present paper be cited in any published work using the 
TRIQS library directly (e.g. for data analysis) or indirectly (e.g. through TRIQS based 
applications). In the latter case, this citation should be added to the citations already 
requested by the application. This helps the TRIQS developers to better keep track of 
projects using the library and provides them guidance for future developments. 

5. Programming Requirements 

TRIQS can be used at different levels of expertise, starting from basic Python inter¬ 
active usage to development of cutting-edge mixed Python/C++ high-performance and 
massively parallel codes, and in pure C++. 

Most objects, in particular the Green’s functions, have a rich Python interface, al¬ 
lowing one to easily plot and manipulate them. For example, simple operations such 
as value assignment, inversion or output to and input from HDF5 files are all one-line 
operations, as shown in the examples below. 

At the C++ level, the required knowledge to make efficient use of the library is min¬ 
imised. The systematic usage of modern C++ (C++11 and C++14) very often lead to 
simpler syntax than old C++. The library often favours a “functional style” programming 
and the simplest possible constructions for the C++ user. To fully exploit the capabilities 
of the library, some understanding of the basic notions of generic programming, such as 
concepts and templates, is helpful, but not required. More traditional object-oriented 
notions of C++ such as inheritance or dynamical polymorphism (virtual functions) are 
not necessary to use the TRIQS library. 


6. Portability 

TRIQS is written in modern C++, i.e. using the C++11 ISO standard. The motivation 
for this choice is twofold: first, we encourage the users of the library to benefit from 
the new features of C++, in particular those which produce much simpler code (e.g. auto, 
for auto loop or lambdas). Second, it dramatically reduces the cost of implementing 
and maintaining the library itself, since many of the new C++ features are designed to 
facilitate the use of the metaprogramming techniques needed to implement high-level, 
high-performance libraries. 

As a result, TRIQS requires a C++11 standard-compliant compiler. The documenta¬ 
tion provides an updated list of tested compilers. When it is available, we recommend 
using a C++14 compiler for development, in particular to get simpler error messages. 

At the Python level, we use the 2.7 versions of Python. Support for Python 3 is 
planned for later releases. We use the binary hierarchical data format (HDF5) to guar¬ 
antee portability of user generated data in binary form. 
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7. TRIQS in two examples 

Here we illustrate the use of the library for two typical tasks encountered in many- 
body physics; this should give a flavour of the possibilities offered by the library. The first 
example is a complete DMFT computation implemented in Python, using a continuous- 
time quantum Monte Carlo solution of the impurity model (which is presented in Ap¬ 
pendix A). The second example illustrates the manipulation of Green’s functions in C++ 
within TRIQS. 

7.1. A DMFT computation in one page of IPython 

This example requires the CT-INT tutorial application to be installed. The installa¬ 
tion procedure is described in Sec. 9.2. 

Fig. 2 shows a screenshot of an IPython notebook implementing a DMFT self- 
consistency loop. The essential steps are to load the solver module, set parameters 
and an initial guess for the Green’s function and to loop over the DMFT iterations. 
The solver module, for which performance is critical, is written in C++ but used from 
Python. This notebook is available in the sources of the ctint_tutorial application (in 
the examples subdirectory) along with the corresponding python script that is suitable 
for parallel execution. 

The Python framework is highly flexible. In this example, we exchange the random 
number generator in the final DMFT iteration to consolidate the result. We could also 
turn on additional measurements, or implement more sophisticated stopping criteria for 
the loop. Using the IPython notebook, results can be plotted and analysed interactively. 
As a minimal example, we load the HDF5 archive from the disk and plot the imaginary 
part of the Green’s functions in the second cell of the notebook. Parameters used in 
the calculation can easily be saved to file and retrieved for data analysis. On a parallel 
machine, the first part of the script is executed in Python (without the notebook), on 
multiple cores, while the analysis can still be done in the IPython notebook. 

Within this framework, DMFT can readily be explored and practised by non-experts. 
The major part of the calculation is of course performed by the solver module. In 
this example, we have used the interaction-expansion continuous-time quantum Monte 
Carlo (CT-INT) solver. We can easily switch to a more sophisticated hybridisation- 
expansion continuous-time quantum Monte Carlo algorithm (CT-HYB) solver by loading 
the appropriate solver module instead. The complete listing of the C++ implementation of 
the solver module is given and explained in Appendix A, as a more detailed illustration 
of the library’s features. 

7.2. Easy manipulation of Green’s functions in C++ 

Here, we show how to compute in C++ a hybridisation function A(r) in imaginary 
time given the bare dispersion of a two-dimensional square lattice with nearest neighbour 
hopping, at chemical potential p. This is a typical task which is usually performed at 


6 


Figure 2: Screenshot of an IPython notebook executing a DMFT loop using the CT-INT 
solver. 
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from pytriqs.gf.local import * 
from pytriqs.archive import * 
import pytriqs.utility.mpi as mpi 

from pytriqs.applications.impurity_solvers.ctinttutorial import CtintSolver 
from pytriqs.plot.mpl_interface import oplot 


# Parameters 
U = 2.5 
mu = U/2.0 
half_bandwidth=l .0 
beta =40.0 
n_iw =128 
n_cycles = 10000 
delta =0.1 
n_iterations =21 


# Hubbard interaction 

# Chemical potential 

# Half bandwidth (energy unit) 

# Inverse temperature 

# Number of Matsubara frequencies 

# Number of MC cycles 

# delta parameter 

# Number of DMFT iterations 


S = CtintSolver(beta, n_iw) # initialize the solver 

S.G_iw « Semicircular(half_bandwidth) # Initialize the Green's function 

with HDFArchive(" dmft_bethe.output. h5", 'w' ) as A: 

A[' n_iterations'] = n_iterations # Save a parameter 

for it in range (n_iterations): # DMFT loop 
for name, GO in S.G0_iws 

GO «= inverse(iOmegan + mu - (half_bandwidth/2 . 0) **2 * S.G_iw[name] ) / Set GO 
# Change random number generator on final iteration 

randomname = 'mtl9937' if it<n_iterations-l else ' lagged_fibonaccil9937' 

S.solve(U, delta, n_cycles, random_name=random_name) # Solve the impurity problem 

G_sym = (S.G_iw[ ’up' ]+S.G_iw[ 'down' ])/2 # Impose paramagnetic solution 
S.G_iw « G_sym 

if mpi.is_master_node(): 

A['G%i'%it] = G_sym # Save G from every iteration to file 

Starting on 1 Nodes at : 2014-10-22 18:06:12.886058 


A = HDFArchive(" dmftbethe.output.h5" ,' r' ) # Open file 
n_iterations = A[' niterations' ] # Load a parameter 
for it in range (niterations): # Use parameter in the analysis 
if not it%2: oplot(A[ 'G%i' %it].imag, '-o', name='Im G%i'%it) # Plot every second result 
xlim(0,5) 


Out[2] : (0, 5) 
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the beginning of a DMFT calculation. The necessary steps are the following: 
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The sum over k = ( k x ,k y ) is taken over the Brillouin zone, uj n is a fermionic Matsub- 
ara frequency and /? is the inverse temperature. Using the library, these equations are 
implemented as follows: 


Listing 1 Computing the hybridisation 


#include <triqs/gfs.hpp> 1 

using namespace triqs::gfs; 2 

using namespace triqs::lattice; 3 

int mainO { 4 

5 

double beta = 10, mu = 0; 6 

int n_f req = 100, n_pts = 100; 7 

8 

// Green’s function on Matsubara frequencies, lxl matrix-valued. 9 

auto Delta_iw = gf<imfreq>{{beta, Fermion, n_freq}, {1, 1}}; 10 

auto Gloc = gf<imfreq>{{beta, Fermion, n_freq}, {1, 1}}; 11 

12 

// Green’s function in imaginary time, lxl matrix-valued. 13 

auto Delta_tau = gf<imtime>{{beta, Fermion, 2 * n_freq + 1}, {1, 1}}; 14 

15 

auto bz = brillouin_zone{bravais_lattice{{ {1, 0}, {0, 1} }}}; 16 

auto bz_mesh = regular_bz_mesh{bz, n_pts}; 17 

18 

triqs :: clef :: placeholder <1 > k_ ; 19 

triqs :: clef :: placeholder <2> iw_ ; 20 

21 

// The actual equations 22 

Gloc(iw_) << sum(l/(iw_ + mu - 2*(cos(k_ [0]) + cos(k_ [1]))) , k_ = bz_mesh) 23 

/bz_mesh.size () ; // (1) 24 

Delta_iw(iw_) << iw_ + mu - 1 / Gloc(iw_); // (2) 25 

Delta_tau() = inverse.fourier(Delta_iw); // (3) 26 

27 

// Write the hybridization to an HDF5 archive 28 

auto file = triqs: :h5 : :file( "Delta.h5" , H5F_ACC_TRUNC) ; 29 

h5_write (file , "Delta_tau", Delta_tau); 30 

h5_write(file, "Delta_iw", Delta_iw); 31 

} 32 


In the implementation of equations (1) and (2) (lines 23-25), we use a compact syntax 
for the assignment to the Green’s function container provided by the TRIQS library (the 
CLEF library, Sec. 8.5). By definition, this is equivalent to assigning the evaluation 
of the expression on the right-hand side to the data points of the “Green’s function” 1 
Deita_iw on the left and for all Matsubara frequencies in its mesh. iw_ is a placeholder , 
i.e. a dummy variable standing for all points in the Green’s function’s mesh. 

In line 23, the formal expression made of iw_ and k_ is summed over the values 
of k_, and assigned to the function for each iw_. Note that no copy is made by this 
statement, the computation is inlined by the compiler, as if it was written manually. 
This technique is more concise than writing a for-loop on each variable, reduces the risk 
of errors and simultaneously increases the readability. Moreover, such techniques come 
with no performance penalty (as indicated by our tests on several standard compilers). 

In line, we assign the inverse Fourier transform of A(w) to A(r) [Eq. (3)]. Note 
that the high-frequency expansion is part of the Green’s function container and so is 
automatically computed in lines 23 and 25 (see Sec. 8.2 for details). It is used to properly 
treat the discontinuity in the Fourier transformation in line 26. 

Finally, the Green’s functions are stored in an HDF5 file with a simple interface, in a 
portable manner. The storage conventions are detailed in the reference documentation. 


1 Here we refer the Green’s function containers and objects representing hybridisation functions or 
functions with the same signature simply as “Green’s functions”. 
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This example is interesting for two reasons: firstly, it shows that the TRIQS library 
performs a lot of low-lying operations. There is no need to reimplement them and the 
user can concentrate on the physics; secondly, it shows that one can write quite complex 
operations concisely, which is necessary in order to write readable codes. 


8. Library components 

In this section, we provide an overview of the TRIQS library components. We il¬ 
lustrate them either with small examples or in the CT-INT impurity solver example 
presented in Appendix A. The description is neither meant to be complete nor ex¬ 
haustive; the online reference documentation (http://ipht.cea.fr/triqs) will fill the 
gaps. 

8.1. Multidimensional arrays (C++) 

TRIQS provides its own multidimensional arrays, with an emphasis on flexibility, 
performance and the Python interface. It is a fundamental building block for higher- 
level containers, such as Green’s functions. Listing 2 below illustrates some of their 
features. 


Listing 2 Array / matrix example 

#include <triqs/arrays.hpp> 
using namespace triqs::arrays; 
int main () { 

auto a = matrix <double >(2, 2); // Declare a 2x2 matrix of double 

auto b = array< double , 3>(5, 2, 2); // Declare a 5x2x2 array of double 

auto c = array< double , 2> {{1,2,3}, {4,5,6}}; // 2x3 array, with initialization 

triqs::clef::placeholder<0> i_; 
triqs::clef::placeholder<1> j_; 
triqs::clef ::placeholder<2> k_ ; 

// Assignment of values using CLEF 

a(i_, j.) << i. + j_; 

b(i_, j_, k_) << i_ * a(k_, j_); 

std::cout << "a = " << a << std::endl; // Printing 

matrix <double > i = inverse(a); // Inverse using LAPACK 

double d = determinant(a); // Determinant using LAPACK 

auto ac = a; // Make a copy (the container is a regular type) 

ac=a*a+2* ac; // Basic operations (uses BLAS for matrix product) 

b(0, range(), range()) = ac; // Assign ac into partial view of b 

// Writing the array into an HDF5 file. 

auto f = triqs::h5::file ("a_file.h5" , H5F_ACC_TRUNC); 

h5_write(f, "a", a); 

auto m = max_element(abs(b)); // maximum of the absolute value of the array. 

// A more "functional" example: compute the norm sum_{i,j} I A _{ij}I 
auto lambda = [] (double r, double x) { return r + std::abs(x); }; 

auto norm = fold(lambda)(a, 0); 

} 


The library provides three types of containers: array (for multidimensional arrays), 
matrix and vector with the following main characteristics: 
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• Regular-type semantics: Just like std::vector, these containers have regular- 
type semantics. 

• Views: Each container has a corresponding view type (e.g. array_view) to e.g. work 
on slices and partial views. 

• CLEF: The containers are compatible with CLEF (Sec. 8.5 and Listing 2) for fast 
assignment techniques. 

• Python interface: These containers can be easily converted to and from Python 
NumPy arrays. 

• Interface to HDF5: See Sec. 8.6 and Listing 2. 

• Arithmetics: Arithmetic operations are implemented using expression templates 
for optimal performance. 

• BLAS/LAPACK: A BLAS/LAPACK interface for matrices and vectors is pro¬ 
vided for the most common operations. 

• STL compatible iterators: The containers and views can be traversed using 
such iterators, or with simple foreach constructs. 

• Optionally, a (slower) debug mode checks for out-of-range operations. 


8.2. Green’s functions (C++ and Python) 

The library provides a special set of containers that allow one to store and manip¬ 
ulate the various Green’s functions used in the quantum many-body problem and its 
algorithms. They are defined on meshes for various domains, they are tensor-, matrix-, 
or scalar-valued and can be block-diagonal. 

Domains currently implemented include real and imaginary frequencies, real and 
imaginary times, Legendre polynomials, and Brillouin zones. Multiple variable Green’s 
functions are also part of the library, but are restricted to C++14 mode only and are of 
alpha quality in release 1.2. We will not discuss them in this paper. 

Green’s functions optionally include a description of their high-frequency behaviour 
in terms of their moments. Storing this information is important for several operations 
(e.g. Fourier transformation, frequency summation) where the high-frequency behaviour 
needs to be treated explicitly. Being part of the object, the singularity is consistently 
recomputed in all arithmetic operations so that the user need not work out the high- 
frequency asymptotics. Listing 3 illustrates some basic usage of Green’s functions, while 
a Python example has been given above (Fig. 2). 


Listing 3 Green’s function example 

#include <triqs/gfs.hpp> 
using namespace triqs; 
using namespace triqs::gfs; 
using namespace triqs::lattice; 
int main () { 
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double beta = 10; 
int n_freq = 1000; 

clef::placeholder<0> iw_ ; 
clef::placeholder<1> k_ ; 

// Construction of a 1x1 matrix-valued fermionic gf on Matsubara frequencies, 
auto g_iw = gf<imfreq>{{beta, Fermion, n_freq}, {1, 1}}; 

// Automatic placeholder evaluation 
g_iw(iw_) << 1 / (iw_ + 2); 

// Inverse Fourier transform to imaginary time 

auto g_tau = gf<imtime>{{beta , Fermion, 2 * n_freq + 1}, {1, 1}}; 

g_tau() = inverse_fourier(g_iw); // Fills a full view of g_tau with FFT result 

// Create a block Green’s function composed of three blocks, 

// labeled a,b,c and made of copies of the g_iw functions, 
auto G_iw = make_block_gf({ "a" , "b", "c"}, {g_iw, g_iw, g_iw}); 

// A multivariable gf: G(k,omega) 

auto bz = brillouin_zone{bravais_lattice{{{l, 0}, {0, 1}}}}; 

auto g_k_iw = gf<cartesian_product<brillouin_zone, imfreq>>{ 

{{bz , 100}, {beta, Fermion, n_freq}}, {1, 1}}; 

g_k_iw(k_, iw_) << 1 / (iw_ - 2 * (cos(k_(0)) + cos(k_(l))) - 1 / (iw_ + 2)); 

// Writing the Green’s functions into an HDF5 file. 

auto f = h.5 : : file ("file_g_k_iw .h5" , H5F_ ACC_TRUNC ) ; 

h5_write(f, "g_k_iw", g_k_iw); 

h5_write(f, "g_iw", g_iw); 

h5_write(f, "g_tau", g_tau); 

h5_write(f, "block_gf", G_iw); 

} 


The main characteristics of Green’s functions are: 

• Arithmetics: Just like arrays, Green’s functions implement arithmetic operations 
using expression templates. 

• Quick assignment: The class uses the CLEF component of the TRIQS library 
for quick assignment (see Sec. 8.5 and Listing 3). 

• Python interface: The Green’s functions are easily shared between Python and 
C++, see Sec. 8.8, and can thus be used in conjunction with the Python visualisation 
tools. 

• Fourier transforms: TRIQS provides a simple interface to fast Fourier trans¬ 
forms (FFTW). For Green’s functions the information about the high-frequency 
behaviour is used to avoid numerical instabilities. 

• Interface to HDF5. 


8.3. Monte Carlo tools (C++) 

The TRIQS library provides several classes for writing Metropolis-like (quantum) 
Monte Carlo algorithms. In addition to some basic analysis tools, like binning or jack¬ 
knife, the library mainly contains the mc_generic class that implements the Metropolis 
algorithm (choose a move, try the move, compute Metropolis ratio, reject or accept, 
etc.) in terms of completely generic moves (configuration updates) and measurements. 
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In practice, one just needs to implement the moves and measurements. The only 
requirement is that they must model their respective concepts 2 . For example, the concept 
of a move is given by Listing 4. Note in particular that they do not require inheritance 
or virtual functions, which makes them particularly simple to use. 


Listing 4 Concept of a Monte Carlo move 

struct my_monte_carlo_move { 

// propose a change in the configuration and return the Metropolis ratio 
double attempt(); 

// the move has been accepted: modify configuration 
double accept (); 

// the move has been rejected: undo configuration changes 
void reject () ; 

> ; 

A concrete usage of the class is shown in the CT-INT solver example (Appendix 
A) . The class is particularly convenient for complex Monte Carlo algorithms with several 
moves: the moves are isolated from the implementation of the Metropolis algorithm itself 
and each move can be implemented independently. 

The Monte Carlo is of course automatically MPI-enabled. Furthermore, random num¬ 
ber generators can easily be changed dynamically to ensure there is no subtle correlation 
effect. 

Listing 5 illustrates a basic application of the tools for statistical analysis on a corre¬ 
lated random series. Let us assume we have two long vectors vi and V2 storing (possibly 
correlated) samples of the random variables X and Y and that we wish to compute es¬ 
timates of (X) and (X)/(Y), together with the corresponding error bars. In both cases, 
the correlation between samples has to be removed using a binning procedure. This 
being done, the first computation is quite straightforward, while the second one further 
requires a jackknife procedure to remove the bias introduced by the nonlinearity. In 
TRIQS, all these operations are performed by the following code snippet, using a little 
library similar to e.g. ALPS/alea [7]: 


Listing 5 Statistics: error analysis 

//fill observable with the series 
observable <double > X, Y; 

for(auto const & x : VI) X << x; //VI: a vector of statistical samples 
for(auto const & y : V2) Y << y; //V2: a vector of statistical samples 

std: :cout << "<X> is approximately " << average_and_error(X) << std : :endl; 
std::cout << "<X>/<Y> is approx. " << average_and_error(X/Y) << std::endl; 

x«x fills the observable x (a stack of the samples) with the values x of vi. average_and_error (x) 
computes an estimate of (X) and of the error A(X), while average_and_error(x/Y) com¬ 
putes an estimate of (X)/(Y) and of A (X)/(Y). 

8-4■ Determinant manipulations (C++) 

The manipulation of determinants is central to many Monte Carlo approaches to 
fermionic problems, see e.g. [10, 12, 13, 11]. Several cases can be abstracted from 


2 In the sense of CH—h concepts. 
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the following mathematical problem. Let us consider a function F(x,y ) taking real or 
complex values (the type of the arguments x and y is arbitrary) and the square matrix 
M defined by 

Mij = F(xi,yj), (4) 

for two sets of parameters {aq} and {yj} of equal length. The problem consists in quickly 
updating M and its inverse A following successive insertions and removals of one or 
two lines (labelled by Xi ) and columns (labelled by yj ) using the Sherman-Morrison and 
Woodbury formulas [17, 18]. 

This generic algorithm is implemented in the TRIQS det_manip class, using BLAS 
Level 2 [19, 20] internally. The class provides a simple API, in order to make these 
manipulations as straightforward and efficient as possible. 

For optimal efficiency within a Monte Carlo framework, the modifications to the 
matrices can be done in two steps: a first step which only returns the determinant ratio 
between the matrix before (M) and after the modification (M'), i.e. £ = det M'/ det M 
(which is generally used in the acceptance rate of a Metropolis move) and a second step 
which updates the matrix and its inverse. This computationally more expensive step 
is usually done only if the Monte Carlo move is accepted. An example of this class is 
employed is the CT-INT solver discussed in Appendix A. 

8.5. CLEF (C++) 

CLEF (Compile-time Lazy Expressions and Functions) is a component of TRIQS 
which allows one to write expressions with placeholders and functions, and to write 
quick assignments. For example, the following - quite involved - equation 

- 9° v ° (5) 

can be coded as quickly as (variables with underscores denote placeholders) 

chiO(s_, sp_)(nu_, nup_, om_) << 

beta * (g[s_](nu_) * g[sp_](nup_) * kronecker(om_)) 

- beta * (g[s_](nu_) * g[s_](nu_ + om_) 

* kronecker(nu_, nup_) * kronecker(s_, sp_)); 

This writing is clearly much simpler and less error-prone than a more conventional five¬ 
fold nested for-loop. At the same time, these expressions are inlined and optimised 
by the compiler, as if the code were written manually. The library also automatically 
optimises the memory traversal (the order of for loops) for performance based on the 
actual memory layout of the container chio. 

The CLEF expressions are very similar to C++ lambdas, except that their variables are 
found by name (the placeholder) instead of a positional argument (in calling a lambda). 
This is much more convenient for complex codes. 

The precise definition of the automatic assignment is as follows. Any code of the form 
(e.g. with three placeholders): 

A(i_,j_,x_) << expression; 

where expression is an expression involving placeholders 3 is rewritten by the compiler 
as follows: 


3 For a precise list of what is allowed in expressions, the reader is referred to the reference documen¬ 
tation. 
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triqs_auto_assign ( 

A, [](auto& i, autofe j,auto& x) { 

return eval(expression, i_=i, j_ = j, x_=x); 

> 

); 

where triqs_auto_assign is a free function defined by the container A, which fills the con¬ 
tainer with the result of the evaluation of the lambda, and eval evaluates the expression 
(eval is a function and is part of CLEF). The precise details of this operation, such as 
the memory traversal order, are encoded in this function. The CLEF quick assignment 
mechanism can therefore easily and efficiently be extended to any object of the library. 
The library provides adaptors to allow standard mathematical functions such as cos or 
abs and std: : vector to be used in expressions. User-defined functions and class methods 
can conveniently be made compatible with the CLEF quick assignment through macros. 

8.6. HDF5 (C++ and Python) 

HDF5 is a standard, portable and compact file format, see http://www. hdf group. 
org. Almost all objects in the TRIQS library (including arrays or Green’s functions) can 
be stored in and retrieved from HDF5 files, from C++ and/or Python, with a simple and 
uniform interface. For example, in C++: 

auto a = array< double , 2> {{1,2,3}, {4,5,6}}; 

{ 

auto f = h5: :file(" data.h5" , ’w’); 

h5_write(f, ’a’, a); 

} 

or, the corresponding code in Python: 

a = numpy.array([[1,2,3],[4,5,6]]) 
with HDFArchive( "data.h5" , ’w’) as f: 

f [ ’ a ’ ] = a 

In Python, the HDFArchive behaves in a similarly as a diet. Therefore, one can reload 
a complex object (e.g. a block Green’s function) in a single command in a script. An 
example can be seen in Listing 1. 

An HDF5 file can be seen as a tree whose leaves are “basic” objects (multidimensional 
rectangular arrays, double, integer, strings, ...). More complex objects are usually 
decomposed by the library into a subtree of smaller objects, which are stored in an 
HDF5 subgroup. For example, a block diagonal Green’s function (of type BlockGf) is 
stored with subgroups containing the Green’s functions it is made of; a Green’s function 
is stored as a subgroup containing the array of data, the mesh, and possibly the high 
frequency singularity. This format , i.e. the precise conventions for the names and types 
of the small objects and the storage order of the data in the arrays, is described in the 
reference documentation. The HDF5 files can be read without the TRIQS library from C, 
C++, Fortran, Python codes, the HDF5 command line tools and with any tool supporting 
this format. This enables publishing data and facilitates sharing them across different 
groups and platforms. The HDF5 format is indeed widely used, e.g. by the ALPS project 

m. 

Note also that the HDF5 files written from C++ or Python have exactly the same 
format. Hence one can straightforwardly load some Green’s functions in Python that 
have been computed and written using a C++ code, or vice-versa. 


// some data 

// open the file 
// write to the file 
// closes the file 


14 


8.1. Second-quantized operators (C++ and Python) 

The theories of strongly-correlated electron systems often use a language of second- 
quantized operators to formulate the problems under consideration. The model Hamil¬ 
tonians as well as the observables of interest are routinely written as polynomials of 
fermionic operators c* and c. 

The TRIQS library implements a C++ template class many_body_operator, which ab¬ 
stracts the notion of a second-quantized operator. The purpose of this class is to make 
expressions for second-quantized operators written in the C++ or Python code as close 
as possible to their analytical counterparts. In order to pursue this goal, the class im¬ 
plements the standard operator algebra. The library stores the expression in normal 
order, so it performs automatically basic simplications, for example when an expression 
vanishes. Any operator can be constructed as a polynomial of the elementary operators 
carrying an arbitrary number of integer/string indices (defined at compile time). The 
coefficients of the polynomials may be real, complex or of a user-defined numeric type in 
advanced use scenarios. 

There is also a Python version of the same class (called Operator), specialised for the 
case with real coefficients and the fermionic operators with two indices (this particular 
choice is made for compatibility with the Green’s function component). Anyone writing 
a TRIQS-based many-body solver may benefit from this class. For example, the user of 
the solver could define a model Hamiltonian in a Python script and subsequently pass it 
to the solver: 

from pytriqs.operators.operators import Operator, n, c_dag , c 

# Spin operators 

Sp = c_dag(" up" , 0)*c("dn" ,0) # S_ + 

Sm = c_dag ( " dn " , 0) * c ( " up " , 0) # S_- 

Sz = 0.5* (n( "up" ,0) - n ( " dn " , 0) ) # S_z 

S2 = Sz*Sz + (Sp*Sm + Sm*Sp)/2 # S~2 

# The Hamiltonian of a half-filled Hubbard atom: four equivalent forms 
U = 1.0 

HI = -U/2*(n(" up" , 0) + n ( " dn " , 0) ) + U*n ( " up " , 0) *n ( " dn " , 0) 

H2 = U*(n( "up" , 0) - 0.5) * ( n ( "dn" , 0) - 0.5) - U/4 

H3 = -2.0*U* Sz * Sz 
H4 = -2.0/3.0*U*S2 

print HI, ’\n■ , H2, ’\n’, H3, ■\n’ , H4 

# All four forms are indeed equivalent 

print (HI-H2) . is_zero() and (H2-H3) . is_zero() and (H3-H4) . is_zero() 


8.8. C++/Python wrapping tool 

The tool that glues together the C++ components to Python is a crucial piece of the 
TRIQS project. Indeed the C++/Python architecture of the project is very demanding 
in this aspect: we need to expose diverse components from C++ to Python. These range 
from simple functions to complex objects with methods, overloaded arithmetic operators, 
the interface to HDF5, and so on. The tool must be very flexible, while being as simple 
as possible to use in the most common cases. 

The TRIQS library proposes such a tool in version 1.2. From a simple Python-written 
description of the classes and functions to expose to Python, it generates the necessary 
C wrapping code to build the Python module. Utilities are also included to actually 
compile and setup the modules with cmake. 

In most cases, the process can be fully automatised, using a second tool based on 
the Clang library, which parses the C++ code using libClang and retrieves the descrip¬ 
tion of the classes and functions along with their documentation. As an example, the 
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automatically produced description files for the CT-INT algorithm is provided in the 
Appendix. 

In more complex cases, some information can be added manually to the class de¬ 
scription, for example the fact that the object forms an algebra over the doubles. In 
such a case, by adding a single line to the description file, the tool automatically gener¬ 
ates all the necessary operators for the algebra structure in Python by calling their C++ 
counterparts. 

As a consequence, this tool also allows the TRIQS user to write C++ code directly 
within the IPython notebook and use it immediately , using a so-called “magic cell com¬ 
mand” , in IPython terminology. This is illustrated in Fig. 3. In this case, the command 

Figure 3: Using C++ directly within the IPython notebook 

In [1]: 

In [2]: 


In [3]: 

In [4]: 

Out[4]: 

"/.•/.triqs extracts the prototype of the C++ heiioO function, writes, compiles and loads the 
Python module to be used in the next cell. TRIQS objects, along with STL containers 
(e.g. vector, tuple), can be used as function arguments or return values. 

Using this feature one can tinker with C++ codes directly inside a Python environment, 
without having to set up a C++ project. It is suitable for debugging, quick testing, 
or executing short C++ code. For longer codes, it is better to set up a Python/C++ 
project along the lines shown for the CT-INT in the Appendix. Note that this feature is 
experimental in release 1.2 and currently limited to a single C++ function per cell (even 
though generalisation is quite straightforward). 

9. Getting started 

Detailed information on installation can be found on the TRIQS website and current 
issues and updates are available on GitHub. 

9.1. Obtaining TRIQS 

The TRIQS source code is available publicly and can be obtained by cloning the 
repository on the GitHub website https://github.com/TRIQS/triqs. As the TRIQS 
project is continuously evolving, we recommend that users always obtain TRIQS from 
GitHub. Fixes to possible issues are also applied to the GitHub source. 


%reload_ext pytriqs.magic 
%%triqs 

#include <string> 

std:sstring hello(){ return "Hello, world!"; } 

print hello() 

Hello, world! 

type (hello()) 
str 
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9.2. Installation 

Installing TRIQS is straightforward. We use the cmake tool to configure, build and 
test the library. Assuming that all dependencies have been installed (refer to the online 
documentation), the library is simply installed by issuing the following commands at the 
shell prompt: 

$ git clone https://github.com/TRIQS/triqs.git src 
$ mkdir build_triqs && cd build_triqs 
$ cmake ../src 
$ make 
$ make test 
$ make install 

By default, the installation directory INSTALL_DIR will be located inside the build di¬ 
rectory. Further installation instructions and help on installing the dependencies can be 
found in the online documentation. 

9.3. Usage 

There are different ways of using TRIQS. In the following, we assume that the location 
of the INSTALL_DIR/bin folder is in the search path. We recommend starting with one 
of the interactive IPython notebook examples provided with this paper (see below). The 
interactive IPython notebook is started using the command 

$ ipytriqs_notebook 

which will open the browser and allow one to open an existing or a new notebook. 
Providing a notebook name as an argument will open the notebook directly. 

The IPython example in Fig. 2 uses the CT-INT solver of Appendix A, which is 
shipped as a separate application. Installing external applications is straightforward. 
The CT-INT application, for example, is installed as follows: 

$ git clone https://github.com/TRIQS/ctint_tutorial.git src_ctint 
$ mkdir build_ctint && cd build_ctint 

$ cmake -DTRIQS_PATH=INSTALL_DIR_ABSOLUTE_PATH ,./src_ctint 
$ make 
$ make test 
$ make install 

where INSTALL_DIR_ABSOLUTE_PATH is the (absolute) path to the TRIQS installation 
directory. The application will be installed into the applications subdirectory in this 
TRIQS installation directory. Assuming that INSTALL_DIR_ABSOLUTE_PATH/bin is in 
the UNIX search path, one can then execute the example notebook in Fig. 2. To this 
end, navigate to the examples directory of the ctint_tutorial application sources and 
issue the following command: 

$ ipytriqs_notebook dmft_bethe.ipynb 
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This will load the notebook inside a browser. Individual cells can be executed by pressing 
[Shift+ENTER] (refer to the IPython notebook documentation). The same directory 
contains a Python script to execute the same DMFT loop from the command line, which 
is another mode to use TRIQS that is better suited for long computations on a parallel 
machine. It can be executed by typing 

$ pytriqs dmft_bethe.py 

or in parallel by running, e.g., 

$ mpirun -np 4 pytriqs dmft_bethe.py 

These commands produce a file dmf t_bethe. output .h5. To plot the Green’s function 
from the final iteration, we can launch ipytriqs and type: 

$ ipytriqs 

In [1]: from pytriqs.archive import * 

In [2]: from pytriqs.gf.local import * 

In [3]: from pytriqs.plot,mpl_interface import oplot, pit 
In [4]: A = HDFArchive("dmft_bethe.output,h5","r") 

In [5]: oplot(A["G20"].imag, "-o", name="Im G20") 

In [6]: pit.show() 

As a starting point for developing an external application, we provide a minimal skele¬ 
ton application called hello_world. It can be installed in the same way as the CT-INT 
solver. The C++ examples of this paper, various IPython notebooks and the hello_world 
are provided in a dedicated GitHub repository https: //github. com/TRIQS/tutorials. 
git. 

10. Contributing 

TRIQS is an open source project and we encourage feedback and contributions from 
the user community to the library and the publication of applications based on it. Issues 
should be reported exclusively via the GitHub web site at https : //github. com/TRIQS/ 
triqs/issues. For contributions, we recommend to use the pull request system on the 
GitHub web site. Before any major contribution, we recommend to coordinate with the 
main TRIQS developers. 

11. Summary 

We have presented the TRIQS library, a Toolbox for Research on Interacting Quan¬ 
tum Systems. This open-source computational physics library provides a framework 
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for the quick development of applications in the field of many-body quantum physics. 
Several applications have been built on this library already. They are available at 
https://github.com/TRIQS and will be described in other publications. 
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Appendix A. A sample application: Interaction expansion continuous-time 
quantum Monte Carlo algorithm 

In this Appendix, we present the implementation of a simple interaction expansion 
continuous-time quantum Monte Carlo algorithm (CT-INT). We have used this solver 
in the DMFT IPython example in Fig. 2. We first briefly recall the formalism of the 
CT-INT algorithm before discussing the code. 

Appendix A.l. Formalism 

We consider the following single-orbital impurity action 


* = -£ 


drdr'd a {T)G 0 f .(r - r , )d (7 (r') + / dr^ int (T), 


o o 


whose interaction term is a slightly modified Hubbard term 

U 
2 


Hint = TT £ (jH ~ aSt ) (% ~ a<4 ) 


(A.l) 


(A.2) 
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with 

= \ + (2<5 Sff - 1)<5. (A.3) 

Here 6 is a free small parameter which reduces the sign problem and S srT is a Kronecker 
symbol. This rewriting of the interaction term results in a shift of the chemical potential 
(absorbed in the bare Green’s function Go): ft = p — The a’s only appear in the 
interaction term. 

The CT-INT algorithm consists in expanding the partition function Z = f T>[d, d]e~ s 
in powers of Hint ■ One obtains: 



x £ < T r(n t (ri) - a Slt )... (n t (r fc ) - o^X^Ot) - a Sl4 -) • • • ( 714 .( 7 *) - n sa )) 0 ■ 

Si.. .Sfc ='f“,4, 
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(A.4) 



where T t is the time ordering operator. In the original CT-INT algorithm proposed in 
Ref. [10], a U = 1 — aW = a and a U = a U = 0. This choice can be shown to eliminate 
the sign problem for the half-filled single band Anderson model. Here we sum over the 
indices to make the formulation slightly more symmetric. It has the advantage that the 
non-interacting Green’s function G does not explicitly depend on cds. Note that in the 
case of all a SilJ being the same, the sum over s, produces a factor 2 k , which cancels 
the 2 k in the denominator. The non-interacting Green’s function has no off-diagonal 
up/down terms, so that the average factorises into product of two correlation functions 
for each spin. Let us furthermore introduce time ordering by replacing the integrals over 
the complete time intervals into a product of time-ordered integrals, 


dn... dr k / JJ(n CT (Ti) - a Si<7 ) ) = k\ dr_ 


dr 2 


dr k x 


o '* -1 a ' 0 o o o 

x (Or( T i) - aSlT )'' ’ ( n t in) - a Sfct )) 0 (( ? R( r i) - a SlJ -) ■ ■ ■ OhC 7 *) - a sa )) 0 . (A.5) 

Using Wick’s theorem and the usual definition for the Green’s function 

G%[t) = - (T T d a <j)d a { 0)) Q , (A.6) 

the averages can be represented by determinants. We hence arrive at 

°° f (_7T\fc 

Z = Z °Y, dn...dT k ^^detZ^detZ^. (A.7) 
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The determinants explicitly read 
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where we have defined the Green’s function Gff as 
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We can sample the partition function (A. 7) by defining a Monte Carlo configuration as 
C := {{ti, Si}, ..., {rfe, Sfc}} and the Monte Carlo weight of a configuration according to 
w(C) = \(—U/2) k det |. The Metropolis acceptance rate for an insertion of a vertex 


is 
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while for a removal, it is 
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The Green’s function can be calculated as 


1 SlnZ 




(3 <5A ct (—t) 


(A.12) 


Carrying out the functional derivative and Fourier transforming yields 

MC 

G°{iu n ) = G° 0 {iu n ) - ^(iu }n )) 2 Y,Y,[Dk}ij 1 e^ {T ^ ) signHOMC). (A.13) 

P C ij 

Separating the Monte Carlo weight, we need to accumulate 

1 MC 

M%iu n ) = x sign[o;(C)], 

P C ij 
MC 

z = 5^sign[o;(C)], 
c 

From M , we can compute the Green’s function as follows [13]: 

G a {iu] n ) = Gg(iw n ) + G'o(iw n )A/°'(*w n )Go(iw n ). (A.16) 


(A.14) 
(A.15) 


Appendix A.2. Implementation 

As an example of an application of the library, we discuss here the complete code 
listing of a fully working, parallelized implementation of the weak-coupling CTQMC 
algorithm described above. How this code can be used in an actual computation is 
illustrated in the DMFT example of Sec. 7.1. Through the use of the various components 
of the library, including gf, mc_toois, det_manip and clef, the full implementation takes 
about 200 lines; it comes with a Python interface. Note that this simple implementation 
can easily be extended: further measurements and moves may be added, or it may be 
generalised to multi-orbital case or to a retarded interaction. 

We divide the code into several listings that we discuss briefly. The purpose is to 
give an illustration of the possibilities of the TRIQS library without entering into all the 
details. We start with the main header file (Listing 6) of the code. It mainly defines and 
provides access to the Green’s functions that are used in the code, in particular in the 
main member function solve. 


Listing 6 CT-INT: the header file 

#include <triqs/gfs.hpp> 

#include <boost/mpi.hpp> 

// -The main class of the solver- 

using namespace triqs::gfs; 
enum spin {up, down}; 

class ctint_solver { 

block_gf<imfreq> gO_iw , gOtilde_iw, g_iw , M_iw; 
block_gf<imtime> gOtilde_tau; 
double double.occ, percent_done_, beta; 
int n_matsubara, n_times.slices; 
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public : 


II Accessors of the class 

block_gf_view<imfreq> G0_iw() { return g0_iw; } 
block_gf_view<imtime> G0_tau() { return gOtilde_tau; } 
block_gf_view<imfreq> G_iw() { return g_iw; } 

ctint_solver( double beta,, int n_iw = 1024, int n_tau = 100001); 

II The method that runs the qmc 
void solve (double U, double delta, 

int n_cycles, int length_cycle = 50, int n_warmup_cycles = 5000, 
std ::string random_name = 
int max_time = -1); 


>; 


Listing 7 defines the Monte Carlo configurations through a simple vector of determi¬ 
nants A.8 (instances of the det_manip class). They contain all the necessary information 
to completely determine a configuration C := {{ti, si}, ..., {r^, Sfc}}. The determinants 
are constructed from a function object gObarvtau, also declared in this listing, that is 
used to fill the elements of the matrix A.8. 


Listing 7 CT-INT: define the configurations 


// -The QMC configuration 

// Argument type of gObar 
struct arg_t { 

double tau; // The imaginary time 

int s; // The auxiliary spin 

> ; 


// The function that appears in the calculation of the determinant 
struct g0bar_tau { 
gf<imtime> const &gt; 
double beta, delta; 
int s ; 

double operator ()(arg_t const &x , arg_t const &y) const { 
if ((x.tau == y.tau)) { // G_\sigma(0~-) - \alpha(\sigma s) 
return 1.0 + gt [0] (0 , 0) - (0.5 + (2 * (s == x.s ? 1 : 0) - 1) * delta); 

> 

auto x_y = x.tau - y.tau; 
bool b = (x_y >= 0); 
if (!b) x_y + = beta; 

double res = gt[closest_mesh_pt(x_y)](0 , 0); 

return (b ? res : -res); // take into account antiperiodicity 

> 

>; 

// The Monte Carlo configuration 
struct configuration { 

// M-matrices for up and down 

std::vector<triqs::det_manip::det_manip<g0bar_tau>> Mmatrices; 

int perturbation_order() const { return Mmatrices[up] .size (); } 

configurat ion(block_gf<imtime> &g0tilde_tau, double beta, double delta) { 

// Initialize the M-matrices. 100 is the initial matrix size 
for (auto spin : {up, down}) 

Mmatrices.emplace_back(g0bar_tau{gOtilde_tau[spin], beta, delta, spin}, 100); 

} 

}; 


Now that the configuration is declared, the next step is to define the Monte Carlo 
moves that are going to act on this configuration. In Listing 8, two moves are im- 
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plemented: the insertion of an interaction vertex at a random imaginary time and the 
removal of a randomly chosen vertex. They are described by classes that must model 
the concept of a Monte Carlo move. In other words they must have the three members 
attempt, accept, reject. The attempt method tries a modification of the configuration and 
returns a Metropolis acceptance ratio (e.g. for the insertion this ratio is given by A. 10). 
The Monte Carlo class will use this ratio to decide wether to accept or reject the proposed 
configuration and then call accept or reject accordingly. Note that for efficiency reasons 
the update of the determinants is done in two steps: in the attempt method only the ratio 
of the new to the old determinant is computed (via try_insert). The actual update of the 
full inverse matrix is performed only if the move is accepted (see the compiete_operation 
call in accept). 


Listing 8 CT-INT: define the moves 

// - QMC move : inserting a vertex - 

struct move_insert { 
configuration *config; 

triqs::mc_tools::random_generator &rng; 
double beta , U; 

double attempt() { // Insert an interaction vertex at time tau with aux spin s 
double tau = rng(beta); 
int s = rng (2) ; 

auto k = config->perturbation_order(); 

auto det_ratio = config->Mmatrices[up].try.insert(k, k, {tau, s}, {tau, s}) * 
config->Mmatrices[down].try_insert(k, k, {tau, s}, {tau, s}); 
return -beta * U / (k + 1) * det_ratio; // The Metropolis ratio 

> 

double accept () { 

for (auto &d : config->Mmatrices) d.complete_operat ion (); // Finish insertion 

return 1.0; 

> 

void rej ect () {} 

>; 

// - QMC move : deleting a vertex - 

struct move.remove { 
configuration *config; 

triqs::mc_tools::random_generator &rng; 
double beta , U; 

double attempt () { 

auto k = config->perturbation_order (); 

if (k <= 0) return 0; // Config is empty, trying to remove makes no sense 
int p = rng(k); // Choose one of the operators for removal 

auto det_ratio = config->Mmatrices[up].try_remove(p, p) * 
config->Mmatrices[down] .try_remove(p, p) ; 
return -k / (beta * U) * det_ratio; // The Metropolis ratio 

> 

double accept () { 

for (auto &d : config->Mmatrices) d.complete_operation (); 
return 1.0; 

> 

void rejectQ {} // Nothing to do 

>; 


The measurement of the Green’s function is shown in Listing 9. It is a simple tran¬ 
scription of Eq. A. 15. Again, the measurements are described by classes that obey the 
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concept of a Monte Carlo measurement: they have a method accumulate which is called 
during the Monte Carlo chain and accumulates data, and a coiiect_resuits method that 
is called at the very end of the calculation. Typically the coiiect_resuits MPI-reduces 
the results from several cores in a parallelized calculation. Note that stdl4: :plus in 
lines 34 and 35 is the C++14 version of std: :plus, which does not require a type, and 
which is provided by TRIQS for backward compatibility to C++11. 


Listing 9 CT-INT: define the measures 

// - QMC measurement - 

struct measure_M { 

configuration const *config; // Pointer to the MC configuration 
block_gf<imfreq> &Mw; // reference to M-matrix 

double beta, Z = 0; 

measure_M(configuration const *config_, block_gf<imfreq> &Mw_ , double beta.) 

: config(config_), Mw(Mw_), beta(beta_) { 

Mw() = 0; 

> 

void accumulate( double sign) { 

Z += sign; 

for (auto spin : {up, down}) { 

// A lambda to measure the M-matrix in frequency- 

auto lambda = [this, spin, sign](arg_t const &x, arg_t const &y, double M) { 
auto coeff = std: :exp (-1_j * M_PI * (x.tau - y.tau) / beta); 
auto fact = coeff * coeff ; 

for (auto const &om : this ->Mw [spin] .mesh()) { 
this ->Mw[spin][om](0, 0) += sign * M * coeff; 
coeff *= fact; 

} 

}; 

foreach(config->Mmatrices[spin], lambda); 

> 

} 

void collect.results(boost::mpi::communicator const &c) { 
boost::mpi::all.reduce (c , Mw , Mw , stdl4 : :plus<> ( )); 
boost::mpi::all.reduce (c , Z, Z, stdl4: :plus<>()); 

Mw = Mw / (-Z * beta); 

} 


The above components are put together in the main solver body shown in Listing 10. 
The first part is the constructor that only defines the dimension of the Green’s functions. 
The second part is the solve method that actually runs the Monte Carlo simulation. It 
first constructs the Fourier transform Gq(t) of the non-interacting Green’s function given 
by the user (it is stored in g0_iw). Then a Monte Carlo simulation is created by adding 
the relevant moves and measures. This is done via the add_move and add_measure methods. 
Note that both the moves and the measurements are constructed with a reference to 
the Monte Carlo configuration config. The simulation is launched with start and final 
results are collected at the end of the simulation with coiiect_resuits. In line 50, we 
finally compute the actual Green’s function through a compact CLEF expression that 
implements (A. 15). 


Listing 10 CT-INT: the main solver body 


24 






1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

1 


II -The main class of the solver - 

ctint.solver::ctint_solver( double beta_, int n_iw, int n_tau) : beta(beta_) { 
g0_iw = 

make_block_gf({" up" , "down"}, gf<imfreq>{{beta, Fermion, n_iw}, {1, 1}}); 

g0tilde_tau - 

make_block_gf({" up" , "down"}, gf<imtime>{{beta, Fermion, n_tau}, {1, 1}}); 

g0tilde_iw = g0_iw; 
g_iw = g0_iw ; 

M_iw = g0_iw; 

} 

// The method that runs the qmc 

void ctint_solver::solve (double U, double delta, int n.cycles, int length_cycle, 

int n_warmup_cycles, std::string random_name, 
int max_time) { 

boost::mpi::communicator world; 
triqs::clef::placeholder<0> spin_; 
triqs::clef::placeholder<1> om_; 

for (auto spin : {up, down}) { // Apply shift to g0_iw and Fourier transform 
gOtilde_iw[spin](om.) << 1.0 / (1.0 / g0_iw[spin](om_) - U / 2); 
gOtilde_tau()[spin] = triqs::gfs::inverse_fourier(gOtilde iw[spin]); 

} 

// Rank - specific variables 

int verbosity = (world.rank() = = 0 ? 3 : 0); 
int random_seed = 34788 + 928374 * world.rank(); 

II Construct a Monte Carlo loop 

triqs::mc_tools::mc_generic<double> CTQMC(n_cycles, length_cycle, 

n_warmup_cycles, random_name, 
random_seed, verbosity); 


II Prepare the configuration 

auto config = configuration{gOtilde_tau, beta, delta}; 

II Register moves and measurements 

CTQMC.add_move(move_insert{&config, CTQMC.rngO, beta, 0}, "insertion"); 
CTQMC.add_move(move_remove{&config, CTQMC.rngO, beta, U}, "removal"); 
CTQMC.add_measure(measure_M{&config, M_iw, beta}, "M measurement"); 

// Run and collect results 

CTQMC.start(1.0, triqs::utility::clock_callback(max_time)); 

CTQMC.collect_results(world); 

II Compute the Green function from Mw 

g_iw [spin_](om_) << gOtilde_iw[spin_](om_) + gOtilde_iw[spin_](om_) * 

M_iw[spin_](om_) * 
g0tilde_iw[spin_](om.); 


The listings above give a complete implementation of the CT-INT algorithm in C++: 
the ctint_soiver class is ready to be used from within other C++ programs. It is however 
convenient to control calculations on the Python level. As discussed above, TRIQS 
provides a tool to easily expose C++ to Python. Starting from a descriptor written 
in Python (shown in Listing 11), it automatically generates a C wrapping code that 
constructs the Python modules for the ctint_soiver. The descriptor basically lists the 
C++ elements that need to be exposed to Python. In most cases, this descriptor can 
be generated automatically by a small analysing tool provided with TRIQS. Here the 
script 11 has been generated automatically using this tool. 


Listing 11 CT-INT: Python wrapper descriptor 

# Generated automatically using the command 
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# c++2py.py ../c++/ctint.hpp -p -m pytriqs.applications.impurity_solvers. 

ctint.tutorial -o ctint_tutorial 
from wrap_generator import * 

# The module 

module = module.(full_name = "pytriqs.applications.impurity_solvers. 
ctint.tutorial" , doc = "") 

# All the triqs C++/Python modules 
module.use_module( ’ gf ’) 

# Add here all includes beyond what is automatically included by the triqs 

modules 

module.add.inclnde (" .,/c++/ctint.hpp" ) 

# Add here anything to add in the C++ code at the start, e.g. namespace 

using 

module.add.preamble( . 

using namespace triqs::gfs; 

II II II ^ 


# The class ctint_solver 
c = class.( 

py.type = "CtintSolver" , # name of the python class 

c_type = " ct int .solver" , # name of the C++ class 

) 


c.add.constructor ("""(double beta_, int n_iw = 1024, int n.tau = 100001)"" 

doc = ...... MM.I ) 


c.add_method( """void solve (double U, double delta, int n.cycles, int 

length.cycle = 50, int n.warmup.cycles = 5000, std ::string random.name 
= int max.time = -1)""", 

doc = """ """) 


c . add.property(name = 
getter 
doc = 


G0_iw" , 

cfunction( "block.gf.view<imfreq> G0_iw ()"), 

.. .. .. .. ^ 


c.add.property(name = 
getter 
doc = 


G0_tau " , 

cfunction (" block_gf_view<imtime> G0_tau ()"), 

.. .. .. .. ^ 


c.add.property(name = 
getter 
doc = 


G_iw " , 

cfunction( "block.gf.view<imfreq> G.iw ()"), 

.. .. .. .. ^ 


module.add.class(c) 


module.generate.code() 


After the code has been compiled and installed a new Python module is available 
in pytriqs.applications.impurity.solvers.ctint.tutorial. The Solver Can then be used as 
illustrated in Fig. 2. As this example shows, C++ members like g0_iw can directly be 
initialised from a Python script and the solve method is also accessible. Controlling 
the solver, or any other C++ code directly from Python makes it very easy to change 
parameters, plot results, build flexible control structures around it, etc., without the 
need to recompile the codes. 
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